#ifndef _MGGS_FUNCTION_
#define _MGGS_FUNCTION_


#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <cmath>
#include "Polynomial.h"

#define PI acos(-1)

using namespace std;

double Bfunc(double x);
double Cfunc(double x);
vector<vector<double> >& BgenerateVector(const int N);
vector<vector<double> >& CgenerateVector(const int N);
vector<vector<bool> >& generateBoolVector(const int N, const int Num = 1);
void plotPolynomial(const string &_file, const Polynomial<double>& p, const double left, const double right, const
                    string& legend ,const double dx = 0.01);
void Bplot(const string &_file);
void Cplot(const string &_file);
void Dplot(const string &_file);
void Eplot(const string &_file);
vector<double>& Chebyshev(const int N);


double Bfunc(double x){
     return 1.0/(1+x*x);
}

double Cfunc(double x){
     return 1.0/(1+25*x*x);
}

vector<vector<double> >& BgenerateVector(const int N){
  vector<vector<double> >* result = new vector<vector<double> >;
  vector<double> vecx(N+1);
  vector<double> vecvalue(N+1);
  for (int i = 0 ; i <= N ; i++ ){
    vecx[i] = -5+10.0*i/N;
    vecvalue[i] = Bfunc(vecx[i]);
  }
  result->push_back(vecx);
  result->push_back(vecvalue);
  return *result;
}

vector<vector<double> >& CgenerateVector(const int N){
  vector<vector<double> >* result = new vector<vector<double> >;
  vector<double> vecx = Chebyshev(N);
  vector<double> vecvalue(N);
  for (int i = 0 ; i < N ; i++ ){
    vecvalue[i] = Cfunc(vecx[i]);
  }
  result->push_back(vecx);
  result->push_back(vecvalue);
  return *result;
}


vector<vector<bool> >& generateBoolVector(const int N, const int Num){
  vector<vector<bool> >* result = new vector<vector<bool> >;
  vector<bool> vecbool[Num+1];
  for (int i = 0 ; i <= Num ; i++){
    vecbool[i].resize(N+1,1);
    result->push_back(vecbool[i]);
  }
  return *result;
}

void plotPolynomial(const string &_file, const Polynomial<double>& p, const double left, const double right, const
                    string &legend, const double dx){
     ofstream os(_file,ios::app);
     const vector<double> coe = p.getCoe();
     const int M = coe.size();
     os << "p=[";
     for(int i = 0 ; i < M ; i++)
       os << coe[M-i-1] << ",";
     os << "];\n";
     os << "x=" << left << ":" << dx << ":" << right << ";\n";
     os << "hold on;\n";
     os << "plot(x,polyval(p,x),'DisplayName','" << legend << "');\n";
     os.close();
}

void Bplot(const string &_file){
     ofstream os(_file,ios::app);
     os << "x=-5:0.1:5;\n";
     os << "hold on;\n";
     os << "plot(x,1./(1+x.^2),'*');legend('exact function');\n";
     os.close();
}


void Cplot(const string &_file){
     ofstream os(_file,ios::app);
     os << "x=-1:0.02:1;\n";
     os << "hold on;\n";
     os << "plot(x,1./(1+25*x.^2),'*');legend('exact function');\n";
     os.close();
}

void Dplot(const string &_file){
     ofstream os(_file,ios::app);
     os << "hold on;\n";
     os << "plot([0,13],[81,81]);legend('81 feet per second');\n";
     os.close();
}

void Eplot(const string &_file){
     ofstream os(_file,ios::app);
     os << "hold on;\n";
     os << "plot([28,43],[0,0]);legend('dying line');\n";
     os.close();
}


vector<double>& Chebyshev(const int N){
     vector<double>* result = new vector<double>;
     for (int i = 1 ; i <= N ; i++){
          result->push_back(cos((2*i-1)*PI/(2*N)));
     }
     return *result;
}



#else
//do nothing
#endif
